{
 "metadata": {
  "name": "",
  "signature": "sha256:b6c351e30d91734e44041c14ce90060c9769838ffe14312e94e02099b1366df2"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as N\n",
      "import pandas as pd\n",
      "\n",
      "import sys\n",
      "sys.path.append('~/PythonLib')\n",
      "from GDALOGRUtilities import OGRWriter\n",
      "\n",
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def arc_distance(lon,lat):\n",
      "    \"\"\"Calculate distance in meters between subsequent points.\"\"\"\n",
      "    \n",
      "    R = 6378.e3 # approximate earth radius\n",
      "    deg2rad = N.pi/180.\n",
      "    \n",
      "    lon = N.asarray(lon)\n",
      "    lat = N.asarray(lat)\n",
      "    d = N.zeros(len(lon),dtype=N.float32)\n",
      "    \n",
      "    dlon = (lon[1:] - lon[0:-1])*deg2rad\n",
      "    dlat = (lat[1:] - lat[0:-1])*deg2rad\n",
      "    latmean = N.mean(lat)*deg2rad\n",
      "    \n",
      "    dx = dlon*N.cos(latmean)*R\n",
      "    dy = dlat*R\n",
      "    \n",
      "    d[0:-1] = N.sqrt(dx**2 + dy**2)\n",
      "    \n",
      "    return d"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def get_reach_index_length(distance,threshold=100.):\n",
      "    \"\"\"Get contiguous river reaches by associating points\n",
      "    separated by a distance less than threshold.\"\"\"\n",
      "    \n",
      "    # Declare an index associated with each point\n",
      "    \n",
      "    reach_index = N.zeros(len(distance),dtype=N.int32)\n",
      "    \n",
      "    # Declare a reach length associated with each point\n",
      "    \n",
      "    reach_length = N.zeros(len(distance),dtype=N.float32)\n",
      "    \n",
      "    # Find all of the locations where the distance exceeds the threshold\n",
      "    \n",
      "    break_idx = N.flatnonzero(distance > threshold)\n",
      "    \n",
      "    nreaches = len(break_idx)\n",
      "    total_reach = N.zeros(nreaches,dtype=N.float32)\n",
      "    npoints = N.zeros(nreaches,dtype=N.int32)\n",
      "    \n",
      "    print('Number of reaches found: %d'%nreaches)\n",
      "    \n",
      "    # Iterate through each reach to set the index and compute the length\n",
      "    \n",
      "    istart = 0\n",
      "    for i in range(nreaches):\n",
      "        iend = break_idx[i]\n",
      "        reach_index[istart:iend+1] = i\n",
      "        reach_length[istart+1:iend+1] = N.cumsum(distance[istart:iend])\n",
      "        npoints[i] = iend - istart + 1\n",
      "        total_reach[i] = N.max(reach_length[istart:iend+1])\n",
      "        \n",
      "        istart = iend + 1\n",
      "        \n",
      "    return reach_index, reach_length, break_idx, total_reach, npoints"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def find_bounding_box(lon,lat,break_idx):\n",
      "    \"\"\"Return the bounding box for each reach.\"\"\"\n",
      "    \n",
      "    lon = N.asarray(lon)\n",
      "    lat = N.asarray(lat)\n",
      "    \n",
      "    nbreak = len(break_idx)\n",
      "    latmin = N.zeros(nbreak,dtype=N.float32)\n",
      "    latmax = N.zeros(nbreak,dtype=N.float32)\n",
      "    lonmin = N.zeros(nbreak,dtype=N.float32)\n",
      "    lonmax = N.zeros(nbreak,dtype=N.float32)\n",
      "    \n",
      "    istart = 0\n",
      "    for i in range(nbreak):\n",
      "        iend = break_idx[i]\n",
      "        lonmin[i] = lon[istart:iend+1].min()\n",
      "        lonmax[i] = lon[istart:iend+1].max()\n",
      "        latmin[i] = lat[istart:iend+1].min()\n",
      "        latmax[i] = lat[istart:iend+1].max()\n",
      "        istart = iend + 1\n",
      "        \n",
      "    return lonmin, lonmax, latmin, latmax"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def get_width_stats(width,break_idx):\n",
      "    \"\"\"Return mean and standard deviation for each reach.\"\"\"\n",
      "    \n",
      "    width = N.asarray(width).astype(N.float32)\n",
      "    \n",
      "    nbreak = len(break_idx)\n",
      "    width_mean = N.zeros(nbreak,dtype=N.float32)\n",
      "    width_std = N.zeros(nbreak,dtype=N.float32)\n",
      "    width_min = N.zeros(nbreak,dtype=N.float32)\n",
      "    width_max = N.zeros(nbreak,dtype=N.float32)\n",
      "    \n",
      "    istart = 0\n",
      "    for i in range(nbreak):\n",
      "        iend = break_idx[i]\n",
      "        width_mean[i] = N.mean(width[istart:iend+1])\n",
      "        width_std[i] = N.std(width[istart:iend+1])\n",
      "        width_min[i] = N.min(width[istart:iend+1])\n",
      "        width_max[i] = N.max(width[istart:iend+1])\n",
      "        istart = iend + 1\n",
      "        \n",
      "    return width_mean, width_std, width_min, width_max"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def write_to_ogr(writer,river_df,reach_df):\n",
      "    \"\"\"Write to a OGR supported driver.\"\"\"\n",
      "    \n",
      "    istart = 0\n",
      "    for i in range(len(reach_df)):\n",
      "        record = reach_df.iloc[i]\n",
      "        iend = int(record['break_idx'])\n",
      "        \n",
      "        lon = N.asarray(river_df['long'][istart:iend+1])\n",
      "        lat = N.asarray(river_df['lat'][istart:iend+1])\n",
      "        \n",
      "        field_record = {'break_idx':int(record['break_idx']),\n",
      "          'npoints':int(record['npoints']),\n",
      "          'reach':float(record['reach']),\n",
      "          'lonmin':float(record['lonmin']),\n",
      "          'lonmax':float(record['lonmax']),\n",
      "          'latmin':float(record['latmin']),\n",
      "          'latmax':float(record['latmax']),\n",
      "          'width_mean':float(record['width_mean']),\n",
      "          'width_std':float(record['width_std']),\n",
      "          'width_min':float(record['width_min']),\n",
      "          'width_max':float(record['width_max']),\n",
      "          }\n",
      "        \n",
      "        if (len(lon) > 2) and (record['npoints']>2) :\n",
      "            writer.add_xy_feature(lon,lat,field_record)\n",
      "        \n",
      "        istart = iend + 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dt = {'width':N.int16,'nchannels':N.int8,'reservoir':N.float32,\n",
      "      'long':N.float32,'lat':N.float32}\n",
      "river_df = pd.read_csv('nAmerica_GRWDL.csv',dtype=dt)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# drop missing data\n",
      "\n",
      "river_df = river_df.dropna()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "distance = arc_distance(river_df['long'],river_df['lat'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "threshold = 1000 # Maximum break before the start of a new reach\n",
      "reach_index, reach_length, break_idx, total_reach, npoints = \\\n",
      "get_reach_index_length(distance,threshold=threshold)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Number of reaches found: 8482\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Add the additional columns to the data frame\n",
      "\n",
      "river_df['reach_index'] = reach_index\n",
      "river_df['reach'] = reach_length"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "river_df.head()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>width</th>\n",
        "      <th>nchannels</th>\n",
        "      <th>reservoir</th>\n",
        "      <th>long</th>\n",
        "      <th>lat</th>\n",
        "      <th>reach_index</th>\n",
        "      <th>reach</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>0</th>\n",
        "      <td> 150</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0</td>\n",
        "      <td>-163.897995</td>\n",
        "      <td> 55.023300</td>\n",
        "      <td> 0</td>\n",
        "      <td>   0.000000</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>1</th>\n",
        "      <td> 150</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0</td>\n",
        "      <td>-163.897995</td>\n",
        "      <td> 55.023102</td>\n",
        "      <td> 0</td>\n",
        "      <td>  22.081333</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>2</th>\n",
        "      <td> 150</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0</td>\n",
        "      <td>-163.897995</td>\n",
        "      <td> 55.022800</td>\n",
        "      <td> 0</td>\n",
        "      <td>  55.627972</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>3</th>\n",
        "      <td> 150</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0</td>\n",
        "      <td>-163.899002</td>\n",
        "      <td> 55.022499</td>\n",
        "      <td> 0</td>\n",
        "      <td> 138.378937</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>4</th>\n",
        "      <td> 192</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0</td>\n",
        "      <td>-163.899002</td>\n",
        "      <td> 55.022301</td>\n",
        "      <td> 0</td>\n",
        "      <td> 160.460266</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "<p>5 rows \u00d7 7 columns</p>\n",
        "</div>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "   width  nchannels  reservoir        long        lat  reach_index       reach\n",
        "0    150          1          0 -163.897995  55.023300            0    0.000000\n",
        "1    150          1          0 -163.897995  55.023102            0   22.081333\n",
        "2    150          1          0 -163.897995  55.022800            0   55.627972\n",
        "3    150          1          0 -163.899002  55.022499            0  138.378937\n",
        "4    192          1          0 -163.899002  55.022301            0  160.460266\n",
        "\n",
        "[5 rows x 7 columns]"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "hist(total_reach/1.e3,bins=N.arange(0,100,0.5),log=True);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEDCAYAAAA2k7/eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFrdJREFUeJzt3X9M3Hcdx/HXbVDjZltb0zLLYViAAbfhqqGakJCd1oZU\nV9TVDTDZDFSjNhjnj8T416iZrGiiNtbE+GNNo4aicYbp7OnqgmnSFmLq4iIotOHijf0w6Q/bEhda\n9vEPBO7gfn7ve/e9732ej6QZ9+Hu8/ncd+37Pvf+/PgGjDFGAICyd5vXHQAAFAcBHwAsQcAHAEsQ\n8AHAEgR8ALAEAR8ALEHABwBLEPABwBIFCfjz8/PatWuXnnvuuUJUDwBwoCAB/1vf+pa6uroKUTUA\nwKGsAn5fX5+qqqrU0tKSUB6JRNTU1KSGhgYNDQ1Jkp5//nmFQiFt27bN/d4CABwLZHOWzunTp/W2\nt71Njz32mF566SVJ0uLiohobG3Xq1ClVV1dr165dGh4e1i9+8QvNz89rcnJSb33rW/Wb3/xGgUCg\n4G8EAJBeRTZPam9vVzQaTSibmJhQfX29amtrJUnd3d0aHR3Vk08+KUk6fvy4tm3bRrAHgBKRVcBP\nZm5uTjU1NSuPg8GgxsfHVx5/6lOfyq9nAABXOQ74+Y7c6+vrdfHixbzqAADb1NXV6cKFC45e63iV\nTnV1tWKx2MrjWCymYDCY9esvXrwoYwx/jNETTzzheR9K5Q/XgmvBtUj/J5+BsuOA39raqpmZGUWj\nUS0sLGhkZESdnZ051TEwMKCxsTGnXQAAa4yNjWlgYCCvOrIK+D09PWpra9P09LRqamp07NgxVVRU\n6OjRo+ro6FAoFFJXV5eam5tzanxgYEDhcNhJvwHAKuFwOO+An1UOf3h4OGn53r17tXfv3rw6APGh\nF4drsYprsYpr4Y6s1uEXpOFAQE888YTC4TD/MwEgg7GxMY2NjenQoUNyGrY9DfgeNQ0AvpVP7OS0\nTACwhKcBn1U6AJAdN1bpkNIBAB8hpQMAyIiADwCWIIcPAD5ADh8ALEMOHwCQEQEfACxBDh8AfIAc\nPgBYhhw+ACAjAj4AWIKADwCWIOADgCVYpQMAPuD7VTr339+m06cj2rhxoxddAADf8e0qnX/+84Ku\nXr3qZRcAwBqeBvzbb9/gZfMAYBUmbQHAEgR8ALAEAR8ALFHhZeMLC//R2bNnVVNT42U3AKDkjY2N\n5b2M3dNlmXfeGdTU1BkCPgBkybfLMgEAxUPABwBLEPABwBIEfACwBAEfACxBwAcASxDwAcASBHwA\nsAQ7bQHAB9hpCwCWYactACAjAj4AWIKADwCWIOADgCUI+ABgCQI+AFiCgA8AliDgA4AlCPgAYAnX\nA/4//vEPff7zn9cjjzyin/70p25XDwBwqGBHK7z55pvq7u7WL3/5y+QNc7QCAOSs4Ecr9PX1qaqq\nSi0tLQnlkUhETU1Namho0NDQ0Er5b3/7W33kIx9Rd3e3o04BANyXVcDv7e1VJBJJKFtcXFR/f78i\nkYgmJyc1PDysqakpSdK+fft08uRJHT9+3P0eAwAcyep45Pb2dkWj0YSyiYkJ1dfXq7a2VpLU3d2t\n0dFR/fvf/9YzzzyjN954Qx/4wAfc7i8AwCHH5+HPzc0l5N6DwaDGx8f1wAMP6IEHHsiqjoWF/+g7\n3/mONm/erHA4rHA47LQ7AFCW3DgHf5njgB8IBPJufMOGzfryl7/MpC0ApLB2MHzo0CHHdTlellld\nXa1YLLbyOBaLKRgMOu4IAKCwHI/wW1tbNTMzo2g0qh07dmhkZETDw8M51cEtDgEgO0W7xWFPT4/+\n/Oc/69KlS9q+fbu+8Y1vqLe3VydPntTjjz+uxcVFHThwQF//+tezb5h1+ACQs3zW4XNPWwDwkXwC\nvuOUjhtI6QBAdoqW0ikERvgAkLuCH60AAPA/UjoA4AOkdADAMqR0AAAZEfABwBLk8AHAB8jhA4Bl\nyOEDADIi4AOAJcjhA4APkMMHAMuQwwcAZETABwBLEPABwBIEfACwBKt0AMAHWKUDAJZhlQ4AICMC\nPgBYgoAPAJYg4AOAJQj4AGAJlmUCgA+wLBMALMOyTABARgR8ALAEAR8ALEHABwBLEPABwBIEfACw\nBAEfACzBxisA8AE2XgGAZdh4BQDIiIAPAJYg4AOAJQj4AGAJzydtpRu67baArl277EU3AMBX8pm0\n9XRZpiTNz1/1ugsAYAVSOgBgCQI+AFiCgA8AliDgA4AlCjJpOzo6queee07Xrl3TgQMHtGfPnkI0\nAwDIQUGXZV69elVf/epX9ZOf/GR9w/9fljk//7IkOV5mBAA2KcpZOn19faqqqlJLS0tCeSQSUVNT\nkxoaGjQ0NJTwuyeffFL9/f2OOgYAcFfWAb+3t1eRSCShbHFxUf39/YpEIpqcnNTw8LCmpqZkjNHX\nvvY17d27Vzt37nS90wCA3GWdw29vb1c0Gk0om5iYUH19vWprayVJ3d3dGh0d1alTp/SnP/1J165d\n04ULF/TZz37WzT4DABzIa9J2bm4u4Sz7YDCo8fFxff/739cXvvCFjK9fWPjPys9jY2MKh8P5dAcA\nyo4bNz5ZllfADwQCeTW+YcNm3bx5XZIcBftNm7ZKkq5du6xNm7bq+vUr2rhxC+fyACgb4XA4IT4e\nOnTIcV15Bfzq6mrFYrGVx7FYTMFgMJ8qc3L9+pU1Pxtdv57fhxAAlKu8An5ra6tmZmYUjUa1Y8cO\njYyMaHh4OOvXx6d0AACpuZHayXqVTk9Pj9ra2jQ9Pa2amhodO3ZMFRUVOnr0qDo6OhQKhdTV1aXm\n5uasG9+wYbOjTqdXsZLqAYByEQ6HNTAwkFcdnp+H72Tj1XK+fkmFpFv//9lICuRcHwD4hW/Pw3ea\n0lnO1y8F91uKD/SZMLkLwI/cSOn4coS/tDooPsgv/5w4wo8P7lL8JO/S8/gWAMBvfDvCL5yKuCWj\n8St3sv8mAADlxpcpncxyS/MAQKmzLqWTOFmbPqWTTdlym+T1AfiFNSmdxMla9+tl0xaAcsYdrwDA\nEmWaw18rfq3+qsQUEQCULqty+KvB2Vm+PtsylmoCKGVFueOV14o9Et+0aasCgQDHNAAoG74J+MW2\nOpF7PWnQ5wMBgN9YksNPJllev0KBwAZJN+PKbiX9dsHKHgDFZFUOP37nbCFz+KnK4vu3dj6BvD+A\nYrEih19KMs0nkO4BUIoI+C6JD/Kr6Z4rSX8PAF7w1U7bUpYpp0/OH4DXLJ60zcXy6ZuVSpzQTf7c\nZJu8ACAfTNoWcdI2XZkxJuUZ/cvvK/73TPICcIpJWwBARiUS8Ct8PKEZf7OV7G3atNWn7xeAX5XI\npO3SDUv8OaGZ7mYr8bn/REuTuEu/5xx+AMVQIgG/XGW685afP+gA+E2JpHQgkeYBUFgsyywhmdI8\nhbwVI7d5BEpbWS3LTLVkcf1NSkpvWabzZZyZy9Zes0It7WTZKOAPZb0sc/WQMtv4eeUSgFJUsgF/\n+ewZOyRb2rk8ocstGAG4o2QDfvmO7FMHdwAopJIN+OUr1+C+dFOWtR8S6U7f5GROAMkQ8EveLS0d\n2Jb4IZHsCOZ0v8v0AQGg/BHwfSn5atp08x6ZPyAAlDt22vpS/PHLa49uTrezF4DNCPi+l+n4BgBY\nUmI7bVd3mcINFdq0aWvBd86ySxcovLLcabv030oljly93kGbbZnX7ScvW3uDluX/5al2Maf7K7H6\nmqUU0saNW+KW0Lq7S5cPEmC9fHbalmhKh1sEFkN8oHb2msKe9Ml9gAF3sUoHa/j/SIdUp46yPwG2\nI+BbKP2xFfkc6eD8w8LNYHz9+pWs9icU4zjqTG3wIYRiKtEcfunlwcsph5/ta5P91UisJ7fXppNs\njsHpCZ7LH2aZThtN9Tw3ZWqDU0qRq7I+LRP+tHbk6sebuzD6Rrkh4Fsn23n6irwC3dr0Sao0i+Ts\naIdCBuPluuPfA8Ef5aBEV+mgcLJdAXWraEcuOGmnkCt4kq1eYsUQygEj/LKX7DhmZ5KPbosxZnB7\n5dD6a+JkBM+oH37DCL/suXf0QvKReDH2TCyvHHJrdL3+mjgZwTPqh9+4PsKfnZ3Vpz/9aT388MNu\nVw3fqEj42b07lyW/N4C7nIyBVr+BZDPqz2UCO5dvEX6cGC8EvnmlVrBlmQ8//LB+9atfpW6YZZm+\nKIv/67H2Juv51peq3mTLMnNZAprNssx8lq1m279USz/jyzIty8xl2WYxlpn6QbkvdS34ssy+vj5V\nVVWppaUloTwSiaipqUkNDQ0aGhpy1AGUssLsul0egSUfrVes/C7X+hL7Wagdw7n3rxx2L6M8ZBXw\ne3t7FYlEEsoWFxfV39+vSCSiyclJDQ8Pa2pqqiCdhFcKcyP11VUw6+/ktZpfv5lzfYn9LNRN4HPv\nHzekR6nIKuC3t7dry5bEI4snJiZUX1+v2tpaVVZWqru7W6Ojo7p8+bI+97nP6cUXX2TUjyJLPV+Q\n3+i6uGsb3MhB51oH+X87OP6bPDc3p5qampXHwWBQ4+Pj2rp1q374wx+60jkgN6lXJOU3ui7u6a1u\nrP7JtQ6+fdjBccB3Y6VE4g1QxvKuD6XOzRU7flf8a1Go+wssfzPgngWF4caNT5Y5DvjV1dWKxWIr\nj2OxmILBYE51bNiwWTdvXv//o7DTrsA3uB3jquJfi0LtG+DbQWGFw2GFw+GVx4cOHXJcl+N1+K2t\nrZqZmVE0GtXCwoJGRkbU2dnpuCMAgMLKKuD39PSora1N09PTqqmp0bFjx1RRUaGjR4+qo6NDoVBI\nXV1dam5uzqnx9fe0RWkq5VRMvhOq/nxv+U6yOpnUzeb5yZ6X6bVOJ6ltm2geGxvTwMBAXnVwHj4b\nr8r+msVvgPK63XR9Sb6hS+t+v/qc3DZ0pas3n3P7c2k3381lqdp1oz6/8O15+IzwUXhejeBzbdf5\nsRHp72DmT+VwPILb74ERfsmVed2+H8u8br+4Zfl828hmNJ+szI8j/GxG6aU+wi9UP3w7wgcAFA8B\nHwAsQQ4f8AWnB7BlmhtIfxxFstemP/xutd5MR0Snej+pylLdgCe+L9keS+3HOQJy+CVX5nX7fizz\nuv3ilrmxYmhtzj1TDj+XtvI5Dtvt+YRMcwzZ9M/psdRuIIcPAPCMp7c4JKWD8lKhwh+0Fp+CKaVN\nY9n0ZSnN48czd+LPC/LqTCI3ztQhpUN6gmvms5SOk9RKMVI6TtM8ma5PKaR0nLw2nzYyPY+UDgAg\nLQI+AFiCgA8AlmAdPlA0pTTJmlr6s3kKuc6jMNcn8f2sbWN1n0LmEz7Xv9b9df2p9xWwDr/kyrxu\n349lXrfvx7LCTto6fa0bk7aZJnKdTNrm2tZq3aknktO9n1zOA4qX6RTV+OcxaQsASIuADwCWIOAD\ngCXYaQuUlex2vCbfEez2pKkb9a2vY3mna6rnr743JzehCWjjxi0O+rlaR7LdxG7szmWnbcmVed2+\nH8u8bt+PZV63n1+ZN5PKzvue66Rtpn5muqELk7YAgLwR8AHAEgR8ALAEAR8ALEHABwBLcJYOgCJa\nfx6N89eWgmR9yvweU9+jd73lc3ruuGMjZ+mUVpnX7fuxzOv2/Vjmdft+LHNehzs3rUlen9OzgViW\nCQBIi4APAJYg4AOAJQj4AGAJAj4AWIKADwCWIOADgCUI+ABgCXbaAkCCioSfC7e7d6nu9DtuE3fp\nstO2pMq8bt+PZV6378cyr9v3Y5nX7a+WJdu5m+1uXnbaAgCyQsAHAEsQ8AHAEgR8ALAEAR8ALEHA\nBwBLEPABwBIEfACwBAEfACyRyx2EszI/P6+DBw/qLW95i8LhsD75yU+63QQAwAHXR/jPPPOMHnnk\nEf3oRz/Ss88+63b1AACHsgr4fX19qqqqUktLS0J5JBJRU1OTGhoaNDQ0JEmam5tTTU2NJOn22293\nubsAAKeyCvi9vb2KRCIJZYuLi+rv71ckEtHk5KSGh4c1NTWlYDCoWCwmSXrzzTfd7zEAwJGsAn57\ne7u2bNmSUDYxMaH6+nrV1taqsrJS3d3dGh0d1UMPPaRf//rXOnjwoDo7OwvSaQBA7hxP2sanbiQp\nGAxqfHxcd9xxh55++mlXOgcAcI/jSdt8bwpQV1cXdxa+tHruc8DHZV6378cyr9v3Y5nX7fuxzOv2\nV39OjJ2py1K9tq6uTk45HuFXV1ev5OolKRaLKRgMZv36CxcuOG0aAOCA4xF+a2urZmZmFI1GtbCw\noJGREXL2AFDCsgr4PT09amtr0/T0tGpqanTs2DFVVFTo6NGj6ujoUCgUUldXl5qbmwvdXwCAU6bI\nTp48aRobG019fb05fPhwsZv31L/+9S8TDodNKBQy9957rzly5IgxxphLly6ZD33oQ6ahocHs2bPH\nXLlyxeOeFs+tW7fMzp07zYMPPmiMsfdaXLlyxezfv980NTWZ5uZmc+7cOWuvxeDgoAmFQua+++4z\nPT095o033rDmWvT29prt27eb++67b6Us3XsfHBw09fX1prGx0fzhD3/IWH9Rz9JJtXbfFpWVlfru\nd7+rv//97zp37px+8IMfaGpqSocPH9aePXs0PT2t3bt36/Dhw153tWiOHDmiUCi0Mmll67X44he/\nqA9/+MOamprS3/72NzU1NVl5LaLRqH784x/r/Pnzeumll7S4uKgTJ05Ycy2S7XlK9d4nJyc1MjKi\nyclJRSIRHTx4MPPep4J8TKVw5swZ09HRsfL4qaeeMk899VQxu1BSPvrRj5rnn3/eNDY2mtdee80Y\nY8yrr75qGhsbPe5ZccRiMbN7927zwgsvrIzwbbwWV69eNXffffe6chuvxaVLl8w999xjLl++bG7e\nvGkefPBB88c//tGqazE7O5swwk/13gcHBxOyJB0dHebs2bNp6y7qCD/Z2v25ublidqFkRKNR/fWv\nf9X73/9+vf7666qqqpIkVVVV6fXXX/e4d8XxpS99Sd/+9rd1222rfw1tvBazs7Patm2bent79d73\nvlef+cxnND8/b+W12Lp1q77yla/oXe96l3bs2KG3v/3t2rNnj5XXYlmq9/7KK68krIzMJp4WNeDn\nu3a/XNy4cUP79+/XkSNHtHHjxoTfBQIBK67T7373O23fvl3vec97ZIxJ+hxbrsWtW7d0/vx5HTx4\nUOfPn9edd965LmVhy7W4ePGivve97ykajeqVV17RjRs39POf/zzhObZci2QyvfdM16WoAT/ftfvl\n4ObNm9q/f78effRRfexjH5O09Kn92muvSZJeffVVbd++3csuFsWZM2f07LPP6u6771ZPT49eeOEF\nPfroo1Zei2AwqGAwqF27dkmSPvGJT+j8+fO66667rLsWf/nLX9TW1qZ3vOMdqqio0EMPPaSzZ89a\neS2Wpfo3sTaevvzyy6qurk5bV1EDvu1r940xOnDggEKhkB5//PGV8s7OTh0/flySdPz48ZUPgnI2\nODioWCym2dlZnThxQh/84Af1s5/9zMprcdddd6mmpkbT09OSpFOnTunee+/Vvn37rLsWTU1NOnfu\nnP773//KGKNTp04pFApZeS2Wpfo30dnZqRMnTmhhYUGzs7OamZnR+973vvSVuT3hkMnvf/97c889\n95i6ujozODhY7OY9dfr0aRMIBMz9999vdu7caXbu3GlOnjxpLl26ZHbv3l32S85SGRsbM/v27TPG\nGGuvxYsvvmhaW1vNu9/9bvPxj3/cXL161dprMTQ0tLIs87HHHjMLCwvWXIvu7m7zzne+01RWVppg\nMGiefvrptO/9m9/8pqmrqzONjY0mEolkrD9gTIoEKgCgrHBPWwCwBAEfACxBwAcASxDwAcASBHwA\nsAQBHwAsQcAHAEsQ8AHAEv8DzuKV8rOfXYAAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10942d310>"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get the bounding box\n",
      "\n",
      "lonmin, lonmax, latmin, latmax = find_bounding_box(river_df['long'],\n",
      "                                                   river_df['lat'],\n",
      "                                                   break_idx)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "width_mean, width_std, width_min, width_max = \\\n",
      "get_width_stats(river_df['width'],break_idx)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Place the reach results in a data frame\n",
      "\n",
      "reach_df = pd.DataFrame({'break_idx':break_idx,\n",
      "                         'npoints':npoints,\n",
      "                         'reach':total_reach,\n",
      "                         'lonmin':lonmin,\n",
      "                         'lonmax':lonmax,\n",
      "                         'latmin':latmin,\n",
      "                         'latmax':latmax,\n",
      "                         'width_mean':width_mean,\n",
      "                         'width_std':width_std,\n",
      "                         'width_min':width_min,\n",
      "                         'width_max':width_max,\n",
      "                         })"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "reach_df.head()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>break_idx</th>\n",
        "      <th>latmax</th>\n",
        "      <th>latmin</th>\n",
        "      <th>lonmax</th>\n",
        "      <th>lonmin</th>\n",
        "      <th>npoints</th>\n",
        "      <th>reach</th>\n",
        "      <th>width_max</th>\n",
        "      <th>width_mean</th>\n",
        "      <th>width_min</th>\n",
        "      <th>width_std</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>0</th>\n",
        "      <td>  661</td>\n",
        "      <td> 55.023300</td>\n",
        "      <td> 54.847401</td>\n",
        "      <td>-163.845993</td>\n",
        "      <td>-163.904007</td>\n",
        "      <td>  662</td>\n",
        "      <td> 26237.072266</td>\n",
        "      <td> 417</td>\n",
        "      <td> 181.191849</td>\n",
        "      <td> 30</td>\n",
        "      <td> 74.898125</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>1</th>\n",
        "      <td> 1097</td>\n",
        "      <td> 54.712601</td>\n",
        "      <td> 54.608898</td>\n",
        "      <td>-164.561005</td>\n",
        "      <td>-164.638000</td>\n",
        "      <td>  436</td>\n",
        "      <td> 17917.480469</td>\n",
        "      <td> 323</td>\n",
        "      <td> 140.463303</td>\n",
        "      <td> 30</td>\n",
        "      <td> 59.323334</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>2</th>\n",
        "      <td> 1529</td>\n",
        "      <td> 55.031898</td>\n",
        "      <td> 54.941799</td>\n",
        "      <td>-163.692993</td>\n",
        "      <td>-163.753006</td>\n",
        "      <td>  432</td>\n",
        "      <td> 17512.716797</td>\n",
        "      <td> 189</td>\n",
        "      <td>  61.976852</td>\n",
        "      <td> 30</td>\n",
        "      <td> 30.376234</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>3</th>\n",
        "      <td> 3438</td>\n",
        "      <td> 59.993599</td>\n",
        "      <td> 59.937599</td>\n",
        "      <td>-163.610001</td>\n",
        "      <td>-164.121002</td>\n",
        "      <td> 1909</td>\n",
        "      <td> 85866.937500</td>\n",
        "      <td> 560</td>\n",
        "      <td> 194.303299</td>\n",
        "      <td> 30</td>\n",
        "      <td> 99.383278</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>4</th>\n",
        "      <td> 3683</td>\n",
        "      <td> 59.940899</td>\n",
        "      <td> 59.928001</td>\n",
        "      <td>-162.000000</td>\n",
        "      <td>-162.098999</td>\n",
        "      <td>  245</td>\n",
        "      <td> 10809.191406</td>\n",
        "      <td> 379</td>\n",
        "      <td> 246.257141</td>\n",
        "      <td> 30</td>\n",
        "      <td> 46.900105</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "<p>5 rows \u00d7 11 columns</p>\n",
        "</div>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "   break_idx     latmax     latmin      lonmax      lonmin  npoints  \\\n",
        "0        661  55.023300  54.847401 -163.845993 -163.904007      662   \n",
        "1       1097  54.712601  54.608898 -164.561005 -164.638000      436   \n",
        "2       1529  55.031898  54.941799 -163.692993 -163.753006      432   \n",
        "3       3438  59.993599  59.937599 -163.610001 -164.121002     1909   \n",
        "4       3683  59.940899  59.928001 -162.000000 -162.098999      245   \n",
        "\n",
        "          reach  width_max  width_mean  width_min  width_std  \n",
        "0  26237.072266        417  181.191849         30  74.898125  \n",
        "1  17917.480469        323  140.463303         30  59.323334  \n",
        "2  17512.716797        189   61.976852         30  30.376234  \n",
        "3  85866.937500        560  194.303299         30  99.383278  \n",
        "4  10809.191406        379  246.257141         30  46.900105  \n",
        "\n",
        "[5 rows x 11 columns]"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Output data\n",
      "\n",
      "store = pd.HDFStore('nAmerica_GRWDL.h5',mode='w',complevel=9)\n",
      "store['river'] = river_df\n",
      "store['reach'] = reach_df\n",
      "store.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "river_df.to_csv('nAmerica_GRWDL_river_topo.csv',mode='w')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "reach_df.to_csv('nAmerica_GRWDL_reach_topo.csv',mode='w')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Define the output fields and the writer\n",
      "\n",
      "driver = 'ESRI Shapefile'\n",
      "output_file = 'nAmerica_GRWDL_river_topo'\n",
      "\n",
      "fields = {'break_idx':'int',\n",
      "          'npoints':'int',\n",
      "          'reach':'float',\n",
      "          'lonmin':'float',\n",
      "          'lonmax':'float',\n",
      "          'latmin':'float',\n",
      "          'latmax':'float',\n",
      "          'width_mean':'float',\n",
      "          'width_std':'float',\n",
      "          'width_min':'float',\n",
      "          'width_max':'float',\n",
      "          }\n",
      "\n",
      "writer = OGRWriter(output_file,fields=fields,geometry='LineString')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "write_to_ogr(writer,river_df,reach_df)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "writer.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(reach_df['total_reach'],reach_df['npoints'],'.')\n",
      "xlim(-5,5)\n",
      "ylim(-5,5)\n",
      "grid()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "KeyError",
       "evalue": "u'no item named total_reach'",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-24-fc335e8a9519>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mreach_df\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'total_reach'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mreach_df\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'npoints'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'.'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mxlim\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mylim\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mgrid\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/er/anaconda/lib/python2.7/site-packages/pandas/core/frame.pyc\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   1656\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getitem_multilevel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1657\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1658\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getitem_column\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1659\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1660\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_getitem_column\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/er/anaconda/lib/python2.7/site-packages/pandas/core/frame.pyc\u001b[0m in \u001b[0;36m_getitem_column\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   1663\u001b[0m         \u001b[0;31m# get column\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1664\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_unique\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1665\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_item_cache\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1666\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1667\u001b[0m         \u001b[0;31m# duplicate columns & possible reduce dimensionaility\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/er/anaconda/lib/python2.7/site-packages/pandas/core/generic.pyc\u001b[0m in \u001b[0;36m_get_item_cache\u001b[0;34m(self, item)\u001b[0m\n\u001b[1;32m   1003\u001b[0m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcache\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1004\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mres\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1005\u001b[0;31m             \u001b[0mvalues\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_data\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1006\u001b[0m             \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_box_item_values\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalues\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1007\u001b[0m             \u001b[0mcache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/er/anaconda/lib/python2.7/site-packages/pandas/core/internals.pyc\u001b[0m in \u001b[0;36mget\u001b[0;34m(self, item)\u001b[0m\n\u001b[1;32m   2871\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_for_nan_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2872\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2873\u001b[0;31m             \u001b[0m_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mblock\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_find_block\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2874\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mblock\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2875\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/er/anaconda/lib/python2.7/site-packages/pandas/core/internals.pyc\u001b[0m in \u001b[0;36m_find_block\u001b[0;34m(self, item)\u001b[0m\n\u001b[1;32m   3183\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3184\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_find_block\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3185\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_check_have\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3186\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mblock\u001b[0m \u001b[0;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mblocks\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3187\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mitem\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mblock\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/er/anaconda/lib/python2.7/site-packages/pandas/core/internals.pyc\u001b[0m in \u001b[0;36m_check_have\u001b[0;34m(self, item)\u001b[0m\n\u001b[1;32m   3190\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_check_have\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3191\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mitem\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3192\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'no item named %s'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mcom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpprint_thing\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3193\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3194\u001b[0m     def reindex_axis(self, new_axis, indexer=None, method=None, axis=0,\n",
        "\u001b[0;31mKeyError\u001b[0m: u'no item named total_reach'"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bad = total_reach == 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "hist(npoints[bad],bins=100);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "npoints[bad].max()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}